﻿

#include "csrMatrix.hpp"

#include <algorithm>
#include <numeric>

#include "sparseMatrix.hpp"
#include "unap.hpp"
#include "utilityBasicFunction.h"
#include "utilityExceptions.h"

#ifdef OMP
#include <omp.h>
#endif

#ifdef UNATPP
#include "unatpp.hpp"
#endif

#if defined(KERNEL_USE_UNROLL_4)
#define SparseDensePlusDot(sum, r, xv, xi, nnz)                         \
  {                                                                     \
    if (nnz > 0) {                                                      \
      label nnz2 = nnz, rem = nnz & 0x3;                                \
      switch (rem) {                                                    \
        case 3:                                                         \
          sum += *xv++ * r[*xi++];                                      \
        case 2:                                                         \
          sum += *xv++ * r[*xi++];                                      \
        case 1:                                                         \
          sum += *xv++ * r[*xi++];                                      \
          nnz2 -= rem;                                                  \
      }                                                                 \
      while (nnz2 > 0) {                                                \
        sum += xv[0] * r[xi[0]] + xv[1] * r[xi[1]] + xv[2] * r[xi[2]] + \
               xv[3] * r[xi[3]];                                        \
        xv += 4;                                                        \
        xi += 4;                                                        \
        nnz2 -= 4;                                                      \
      }                                                                 \
      xv -= nnz;                                                        \
      xi -= nnz;                                                        \
    }                                                                   \
  }

#elif defined(KERNEL_USE_UNROLL_2)
#define SparseDensePlusDot(sum, r, xv, xi, nnz)           \
  {                                                       \
    label __i, __i1, __i2;                                \
    for (__i = 0; __i < nnz - 1; __i += 2) {              \
      __i1 = xi[__i];                                     \
      __i2 = xi[__i + 1];                                 \
      sum += (xv[__i] * r[__i1] + xv[__i + 1] * r[__i2]); \
    }                                                     \
    if (nnz & 0x1) sum += xv[__i] * r[xi[__i]];           \
  }
#else
#define SparseDensePlusDot(sum, r, xv, xi, nnz) \
  {                                             \
    label __i;                                  \
    for (__i = 0; __i < nnz; __i++) {           \
      sum += xv[__i] * r[xi[__i]];              \
    }                                           \
  }
#endif

namespace UNAP {

namespace TOOLS {
template <typename Cmpt>
void sortRow(label* col, Cmpt* val, int n) {
  for (int j = 1; j < n; ++j) {
    label c = col[j];
    Cmpt v = val[j];

    int i = j - 1;

    while (i >= 0 && col[i] > c) {
      col[i + 1] = col[i];
      val[i + 1] = val[i];
      i--;
    }

    col[i + 1] = c;
    val[i + 1] = v;
  }
}
}  // namespace TOOLS

template <typename Cmpt>
label csr<Cmpt>::scanRowSize() {
  std::partial_sum(ptr, ptr + nrows + 1, ptr);
  return ptr[nrows];
}
template <typename Cmpt>
void csr<Cmpt>::setSize(label n, label m) {
  nrows = n;
  ncols = m;
  ptr = new label[nrows + 1];
  for (label i = 0; i <= nrows; ++i) ptr[i] = 0;
}
template <typename Cmpt>
void csr<Cmpt>::setNnz(label n) {
  nnz = n;
  col = new label[nnz];
  val = new Cmpt[nnz];
}
template <typename Cmpt>
void csr<Cmpt>::sortRows() {
  for (label i = 0; i < nrows; ++i) {
    const label beg = ptr[i];
    const label end = ptr[i + 1];
    TOOLS::sortRow(col + beg, val + beg, end - beg);
  }
}

template <typename Cmpt>
csr<Cmpt>* csr<Cmpt>::transpose() {
  const label n = nrows;
  const label m = ncols;
  csr<Cmpt>* T = new csr<Cmpt>();
  T->setSize(m, n);
  for (label j = 0; j < nnz; ++j) ++(T->ptr[col[j] + 1]);

  T->scanRowSize();
  T->setNnz(nnz);
  for (label i = 0; i < n; i++) {
    for (label j = ptr[i]; j < ptr[i + 1]; ++j) {
      const label head = T->ptr[col[j]]++;
      T->col[head] = i;
      T->val[head] = val[j];
    }
  }
  std::rotate(T->ptr, T->ptr + m, T->ptr + m + 1);
  T->ptr[0] = 0;
  return T;
}

template <typename Cmpt>
csr<Cmpt>* csr<Cmpt>::product(const csr<Cmpt>& A, const csr<Cmpt>& B,
                              bool sort) {
  csr<Cmpt>* c = new csr<Cmpt>();
  csr<Cmpt>& C = *c;
  C.setSize(A.nrows, B.ncols);
  C.ptr[0] = 0;

#ifdef OMP
#pragma omp parallel
#endif
  {
    std::vector<label> marker(B.ncols, -1);
#ifdef OMP
#pragma omp for
#endif
    for (label ia = 0; ia < A.nrows; ++ia) {
      label Ccols = 0;
      for (label ja = A.ptr[ia], ea = A.ptr[ia + 1]; ja < ea; ++ja) {
        label ca = A.col[ja];
        for (label jb = B.ptr[ca], eb = B.ptr[ca + 1]; jb < eb; ++jb) {
          label cb = B.col[jb];
          if (marker[cb] != ia) {
            marker[cb] = ia;
            ++Ccols;
          }
        }
      }
      C.ptr[ia + 1] = Ccols;
    }
  }
  C.setNnz(C.scanRowSize());

#ifdef OMP
#pragma omp parallel
#endif
  {
    std::vector<label> marker(B.ncols, -1);
#ifdef OMP
#pragma omp for
#endif
    for (label ia = 0; ia < A.nrows; ++ia) {
      label rowBeg = C.ptr[ia];
      label rowEnd = rowBeg;

      for (label ja = A.ptr[ia], ea = A.ptr[ia + 1]; ja < ea; ++ja) {
        label ca = A.col[ja];
        Cmpt va = A.val[ja];

        for (label jb = B.ptr[ca], eb = B.ptr[ca + 1]; jb < eb; ++jb) {
          label cb = B.col[jb];
          Cmpt vb = B.val[jb];
          if (marker[cb] < rowBeg) {
            marker[cb] = rowEnd;
            C.col[rowEnd] = cb;
            C.val[rowEnd] = va * vb;
            ++rowEnd;
          } else {
            C.val[marker[cb]] += va * vb;
          }
        }
      }
      if (sort) TOOLS::sortRow(C.col + rowBeg, C.val + rowBeg, rowEnd - rowBeg);
    }
  }
  return c;
}

namespace PATCH {
inline static bool compareCount(const label& s) { return s > 0; }
inline static bool compareIdxIJ(const IdxIJ& a, const IdxIJ& b) {
  if (a.i_ != b.i_) return (a.i_ < b.i_);
  return (a.j_ < b.j_);
}
template <typename Cmpt>
inline static bool compareTri(const Triplet<Cmpt>& a, const Triplet<Cmpt>& b) {
  if (a.r_ != b.r_) return (a.r_ < b.r_);
  return (a.c_ < b.c_);
};
class CompareVec {
 private:
  const label* vec_;

 public:
  CompareVec(label* v) : vec_(v) {}
  bool operator()(label i, label j) { return vec_[i] < vec_[j]; }
};

}  // namespace PATCH

template <typename Cmpt>
inline CSRMatrix<Cmpt>::CSRMatrix() : Matrix<Cmpt>(), diagPtr_(NULL) {}

template <typename Cmpt>
inline CSRMatrix<Cmpt>::CSRMatrix(Communicator* otherComm)
    : Matrix<Cmpt>(otherComm), diagPtr_(NULL) {}

template <typename Cmpt>
inline CSRMatrix<Cmpt>::CSRMatrix(label n, label nnz, bool symm)
    : SparseMatrix<Cmpt>(n, nnz, symm), diagPtr_(NULL) {
  this->ptr_ = new label[n + 1];
  this->ind_ = new label[nnz];
  this->val_ = new Cmpt[nnz];
}

template <typename Cmpt>
CSRMatrix<Cmpt>::CSRMatrix(label beginRow, label n, label nnz, label* cooI,
                           label* cooJ, Cmpt* cooV, Communicator* comm,
                           bool symm)
    : Matrix<Cmpt>(comm),
      diagPtr_(NULL),
      beginRow_(beginRow),
      beginCol_(beginRow),
      localRows_(n),
      localCols_(n),
      localnnz_(nnz) {
  const int nProcs = comm->getMySize();
  label* dims = new label[nProcs];
  comm->allGather("getGlobalDim", &n, sizeof(label), dims, sizeof(label));
  comm->allReduce("sumLocalnnz", &lnnz(), &this->nnz(), 1, COMM_LABEL,
                  COMM_SUM);

  std::vector<Triplet<Cmpt> > Atuple;
  Atuple.reserve(nnz);
  for (label i = 0; i < nnz; i++) {
    Atuple.push_back(Triplet<Cmpt>(cooI[i], cooJ[i], cooV[i]));
  }
  vtxdist().resize(nProcs + 1);
  label index = 0;
  comm->finishTask("getGlobalDim");
  for (label i = 0; i < nProcs; ++i) {
    vtxdist(i) = index;
    index += dims[i];
  }
  vtxdist(nProcs) = index;
  this->gRows_ = index;

  this->ptr() = new label[lRows() + 1];
  this->ind() = new label[lnnz()];
  this->val() = new Cmpt[lnnz()];
  label row = -1, count = 0;
  for (label i = 0; i < lnnz(); i++) {
    this->val(i) = Atuple[i].v_;
    this->ind(i) = Atuple[i].c_;
    label newRow = Atuple[i].r_;
    if (newRow != row) {
      this->ptr(count) = i;
      row = newRow;
      count++;
    }
  }
  this->ptr(lRows()) = lnnz();

  delete[] dims;

  splitOffDiag();
  comm->finishTask("sumLocalnnz");
}

template <typename Cmpt>
CSRMatrix<Cmpt>::CSRMatrix(label* rowPtr, label* colInd, Cmpt* values,
                           label* dist, bool symm, Communicator* comm,
                           bool deepCopy)
    : Matrix<Cmpt>(comm), diagPtr_(NULL) {
  const int P = comm->getMySize();
  const int rank = comm->getMyId();
  localRows_ = dist[rank + 1] - dist[rank];
  localCols_ = localRows_;
  localnnz_ = rowPtr[localRows_] - rowPtr[0];
  beginRow_ = dist[rank];
  beginCol_ = beginRow_;
  vtxdist_.resize(P + 1);
  std::copy(dist, dist + P + 1, vtxdist_.data());
  this->gRows_ = dist[P];
  label lnnz(localnnz_);
  reduceSum(&lnnz, comm);
  this->nnz_ = lnnz;
  this->symm_ = symm;

  if (deepCopy) {
    this->ptr_ = new label[localRows_ + 1];
    std::copy(rowPtr, rowPtr + localRows_ + 1, this->ptr_);
    this->ind_ = new label[localnnz_];
    std::copy(colInd, colInd + localnnz_, this->ind_);
    this->val_ = new Cmpt[localnnz_];
    std::copy(values, values + localnnz_, this->val_);
  } else {
    this->ptr_ = rowPtr;
    this->ind_ = colInd;
    this->val_ = values;
  }
  splitOffDiag();
}

template <typename Cmpt>
CSRMatrix<Cmpt>::CSRMatrix(label rows, label* rowPtr, label* colInd,
                           Cmpt* values, bool symm, bool deepCopy)
    : Matrix<Cmpt>(NULL), diagPtr_(NULL) {
  localRows_ = rows;
  localCols_ = localRows_;
  localnnz_ = rowPtr[localRows_] - rowPtr[0];
  beginRow_ = 0;
  this->gRows_ = rows;
  this->nnz_ = localnnz_;
  this->symm_ = symm;

  if (deepCopy) {
    this->ptr_ = new label[localRows_ + 1];
    std::copy(rowPtr, rowPtr + localRows_ + 1, this->ptr_);
    this->ind_ = new label[localnnz_];
    std::copy(colInd, colInd + localnnz_, this->ind_);
    this->val_ = new Cmpt[localnnz_];
    std::copy(values, values + localnnz_, this->val_);
  } else {
    this->ptr_ = rowPtr;
    this->ind_ = colInd;
    this->val_ = values;
  }
  splitOffDiag();
}

template <typename Cmpt>
CSRMatrix<Cmpt>::CSRMatrix(const CSRMatrix<Cmpt>& Aseq, Communicator* comm)
    : diagPtr_(NULL) {
  this->commcator_ = comm;
  const int rank = comm->getMyId();
  const int P = comm->getMySize();
  this->gRows_ = Aseq.tsize();
  this->nnz_ = Aseq.nnz();
  this->symm_ = Aseq.symm();
  vtxdist_.resize(P + 1);

  vtxdist_[0] = 0;

  label avgColNum = this->gRows_ / P;
  label remainder = this->gRows_ % P;
  for (label i = 0; i < P; ++i) {
    if (i < remainder) {
      vtxdist_[i] = i * (avgColNum + 1);
    } else {
      vtxdist_[i] = i * avgColNum + remainder;
    }
  }
  vtxdist_[P] = this->gRows_;

  beginRow_ = vtxdist_[rank];
  beginCol_ = beginRow_;
  label erow = vtxdist_[rank + 1];
  localRows_ = erow - beginRow_;
  localCols_ = localRows_;
  localnnz_ = Aseq.ptr(erow) - Aseq.ptr(beginRow_);
  label i0 = Aseq.ptr(beginRow_);
  label i1 = Aseq.ptr(erow);
  this->ptr_ = new label[localRows_ + 1];
  this->ind_ = new label[localnnz_];
  this->val_ = new Cmpt[localnnz_];

  std::copy(Aseq.ptr() + beginRow_, Aseq.ptr() + erow + 1, this->ptr_);
  std::copy(Aseq.ind() + i0, Aseq.ind() + i1, this->ind_);
  std::copy(Aseq.val() + i0, Aseq.val() + i1, this->val_);

  for (label r = localRows_; r >= 0; r--) this->ptr_[r] -= this->ptr_[0];
  splitOffDiag();
}

template <typename Cmpt>
CSRMatrix<Cmpt>::CSRMatrix(const CSRMatrix<Cmpt>& Acopy)
    : Matrix<Cmpt>(Acopy),
      diagPtr_(NULL),
      vtxdist_(Acopy.vtxdist_),
      vtxdistCol_(Acopy.vtxdistCol_),
      beginRow_(Acopy.beginRow_),
      beginCol_(Acopy.beginCol_),
      localRows_(Acopy.localRows_),
      localCols_(Acopy.localCols_),
      localnnz_(Acopy.localnnz_),
      offDiagNnz_(Acopy.offDiagNnz_),
      offDiagStart_(Acopy.offDiagStart_),
      spmvBufs_(Acopy.spmvBufs_) {
  if (Acopy.diagPtr_ != NULL) {
    diagPtr_ = new Vector<Cmpt>(*Acopy.diagPtr_);
  }
  this->gRows_ = Acopy.gRows_;
  this->gCols_ = Acopy.gCols_;
  this->nnz_ = Acopy.nnz_;
  this->symm_ = Acopy.symm_;
  this->ptr_ = new label[localRows_ + 1];
  this->ind_ = new label[localnnz_];
  this->val_ = new Cmpt[localnnz_];
  memcpy(this->ptr_, Acopy.ptr_, (localRows_ + 1) * sizeof(label));
  memcpy(this->ind_, Acopy.ind_, (localnnz_) * sizeof(label));
  memcpy(this->val_, Acopy.val_, (localnnz_) * sizeof(Cmpt));
  splitOffDiag();
}

class Comp {
 private:
  const label* r_;
  const label* c_;

 public:
  Comp(const label* r, const label* c) : r_(r), c_(c) {}
  bool operator()(label i, label j) {
    if (r_[i] == r_[j]) return c_[i] < c_[j];
    return r_[i] < r_[j];
  }
};

template <typename Cmpt>
void CSRMatrix<Cmpt>::print() const {
  POUT << "nProcs: " << this->commcator_->getMySize() << std::endl;
  POUT << "rank: " << this->commcator_->getMyId() << std::endl;
  POUT << "vtxdist: ";
  if (vtxdist_.size() != 0) {
    for (label i = 0; i <= this->commcator_->getMySize(); i++)
      POUT << vtxdist_[i] << " ";
  } else {
    POUT << " null";
  }
  POUT << std::endl << "beginRow: " << beginRow_ << std::endl;
  POUT << "beginCol: " << beginCol_ << std::endl;
  this->SparseMatrix<Cmpt>::print();
  POUT << ENDL;
  POUT << localRows_ << " " << localCols_ << " " << localnnz_ << ENDL;
  for (label i = 0; i <= localRows_; i++) {
    for (label j = this->ptr_[i]; j < this->ptr_[i + 1]; j++) {
      POUT << i + 1 << " " << 1 + this->ind_[j] << " " << this->val_[j] << ENDL;
    }
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::splitOffDiag() {
  sortRows();
  offDiagStart_.resize(localRows_);
  if (localRows_ == this->gRows_) {
    for (label row = 0; row < localRows_; row++)
      offDiagStart_[row] = this->ptr_[row + 1];
    return;
  }

  offDiagNnz_ = 0;
#ifdef OMP
#pragma omp parallel for
#endif
  for (label row = 0; row < localRows_; row++) {
    label lo = this->ptr_[row];
    label hi = this->ptr_[row + 1];
    label* firstInd = &this->ind_[lo];
    label* lastInd = &this->ind_[hi];
    Cmpt* firstVal = &this->val_[lo];
    Cmpt* lastVal = &this->val_[hi];
    while (1) {
      while ((firstInd != lastInd) && (*firstInd) >= beginCol_ &&
             (*firstInd) < beginCol_ + localCols_) {
        ++firstInd;
        ++firstVal;
      }
      if (firstInd == lastInd) {
        lastVal--;
        lastInd--;
        break;
      }
      lastVal--;
      lastInd--;
      while ((firstInd != lastInd) && !((*lastInd) >= beginCol_ &&
                                        (*lastInd) < beginCol_ + localCols_)) {
        --lastInd;
        --lastVal;
      }
      if (firstInd == lastInd) break;

      std::iter_swap(firstInd++, lastInd);
      std::iter_swap(firstVal++, lastVal);
    }
    offDiagStart_[row] = lo + std::distance(&this->ind_[lo], firstInd);
    offDiagNnz_ += this->ptr_[row + 1] - offDiagStart_[row];
  }

  const label n = localRows_;
  label nlocal;
  label* ii;
  Cmpt* vv;
  for (label r = 0; r < n; r++) {
    nlocal = offDiagStart_[r] - this->ptr_[r];
    ii = this->ind_ + this->ptr_[r];
    vv = this->val_ + this->ptr_[r];
    std::vector<label> index(nlocal);
    PATCH::CompareVec comp(ii);
    UNAP::iota(index.begin(), index.end(), 0);
    std::sort(index.begin(), index.end(), comp);
    std::vector<label> adj(nlocal);
    std::vector<Cmpt> val(nlocal);
    for (label c = 0; c < nlocal; c++) {
      adj[c] = ii[index[c]];
      val[c] = vv[index[c]];
    }
    memcpy(ii, adj.data(), nlocal * sizeof(label));
    memcpy(vv, val.data(), nlocal * sizeof(Cmpt));
  }
  for (label r = 0; r < n; r++) {
    nlocal = this->ptr_[r + 1] - offDiagStart_[r];
    ii = this->ind_ + offDiagStart_[r];
    vv = this->val_ + offDiagStart_[r];
    std::vector<label> index(nlocal);
    PATCH::CompareVec comp(ii);
    UNAP::iota(index.begin(), index.end(), 0);
    std::sort(index.begin(), index.end(), comp);
    std::vector<label> adj(nlocal);
    std::vector<Cmpt> val(nlocal);
    for (label c = 0; c < nlocal; c++) {
      adj[c] = ii[index[c]];
      val[c] = vv[index[c]];
    }
    memcpy(ii, adj.data(), nlocal * sizeof(label));
    memcpy(vv, val.data(), nlocal * sizeof(Cmpt));
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::sortRows() {
  for (label i = 0; i < localRows_; ++i) {
    const label beg = this->ptr_[i];
    const label end = this->ptr_[i + 1];
    TOOLS::sortRow(this->ind_ + beg, this->val_ + beg, end - beg);
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::setupSpMVBuffers() const {
  if (spmvBufs_.initialized_) return;
  spmvBufs_.initialized_ = true;

  label offDiagNnz = 0;
#ifdef OMP
#pragma omp parallel for reduction(+ : offDiagNnz)
#endif
  for (label r = 0; r < localRows_; r++)
    offDiagNnz += this->ptr_[r + 1] - offDiagStart_[r];

  const int P = this->commcator_->getMySize();
  MPI_Comm comm = this->commcator_->getMpiComm();
  std::vector<label> spmvRind;
  offDiagNnz_ = offDiagNnz;
  spmvRind.reserve(offDiagNnz);
  for (label r = 0; r < localRows_; r++)
    for (label j = offDiagStart_[r]; j < this->ptr_[r + 1]; j++)
      spmvRind.push_back(this->ind_[j]);
  std::stable_sort(spmvRind.begin(), spmvRind.end());
  spmvRind.erase(std::unique(spmvRind.begin(), spmvRind.end()), spmvRind.end());

  if (vtxdistCol_.size() == 0) {
    if (beginRow_ == beginCol_ && localRows_ == localCols_) {
      vtxdistCol_ = vtxdist_;
    } else {
      label dims[P];
      vtxdistCol_.resize(P + 1);
      label index = 0;
      this->commcator_->allGather("getGlobalDim", &localCols_, sizeof(label),
                                  dims, sizeof(label));
      this->commcator_->finishTask("getGlobalDim");
      for (label i = 0; i < P; ++i) {
        vtxdistCol_[i] = index;
        index += dims[i];
      }
      vtxdistCol_[P] = index;
    }
  }

  for (ulabel r = 0, d = 0, last = -1; r < spmvRind.size(); r++) {
    while (spmvRind[r] >= vtxdistCol_[d + 1]) ++d;
    if (last < d) {
      last = d;
    }

    spmvBufs_.rind_[spmvRind[r]] = std::make_pair(d, r);
  }

  spmvBufs_.prbuf_.reserve(offDiagNnz);
  for (label r = 0; r < localRows_; r++)
    for (label j = offDiagStart_[r]; j < this->ptr_[r + 1]; j++)
      spmvBufs_.prbuf_.push_back(std::distance(
          spmvRind.begin(),
          std::lower_bound(spmvRind.begin(), spmvRind.end(), this->ind_[j])));

  std::vector<label> rsizes(P), ssizes(P);
  for (label p = 0, j = 0; p < P; p++)
    while (j < (label)spmvRind.size() && spmvRind[j] < vtxdistCol_[p + 1]) {
      j++;
      rsizes[p]++;
    }

  MPI_Alltoall(rsizes.data(), 1, mpiType<label>(), ssizes.data(), 1,
               mpiType<label>(), comm);

  label npr = std::count_if(rsizes.begin(), rsizes.end(), PATCH::compareCount);
  label nps = std::count_if(ssizes.begin(), ssizes.end(), PATCH::compareCount);
  spmvBufs_.sranks_.reserve(nps);
  spmvBufs_.soff_.reserve(nps + 1);
  spmvBufs_.rranks_.reserve(npr);
  spmvBufs_.roffs_.reserve(npr + 1);
  label osetRecv = 0, osetSend = 0;
  for (label p = 0; p < P; p++) {
    if (ssizes[p] > 0) {
      spmvBufs_.sranks_.push_back(p);
      spmvBufs_.soff_.push_back(osetSend);
      osetSend += ssizes[p];
    }
    if (rsizes[p] > 0) {
      spmvBufs_.rranks_.push_back(p);
      spmvBufs_.roffs_.push_back(osetRecv);
      osetRecv += rsizes[p];
    }
  }
  spmvBufs_.soff_.push_back(osetSend);
  spmvBufs_.roffs_.push_back(osetRecv);
  spmvBufs_.sind_.resize(osetSend);
  std::vector<MPI_Request> req(npr + nps);
  for (label p = 0; p < npr; p++)
    MPI_Isend(spmvRind.data() + spmvBufs_.roffs_[p],
              spmvBufs_.roffs_[p + 1] - spmvBufs_.roffs_[p], mpiType<label>(),
              spmvBufs_.rranks_[p], 0, comm, &req[p]);
  for (label p = 0; p < nps; p++)
    MPI_Irecv(spmvBufs_.sind_.data() + spmvBufs_.soff_[p],
              spmvBufs_.soff_[p + 1] - spmvBufs_.soff_[p], mpiType<label>(),
              spmvBufs_.sranks_[p], 0, comm, &req[npr + p]);

  std::vector<label> nnz(localRows_);
  spmvBufs_.compRow_.resize(1 + localRows_);
  spmvBufs_.rowIndex_.resize(localRows_);
  spmvBufs_.compRow_[0] = 0;
  label count = 0;
  for (label i = 0; i < localRows_; i++) {
    label nnzRi = this->ptr_[i + 1] - offDiagStart_[i];
    if (nnzRi != 0) {
      spmvBufs_.compRow_[count + 1] = spmvBufs_.compRow_[count] + nnzRi;
      spmvBufs_.rowIndex_[count] = i;
      count++;
    }
  }
  spmvBufs_.compRow_.resize(count + 1);
  spmvBufs_.rowIndex_.resize(count);

  MPI_Waitall(req.size(), req.data(), MPI_STATUSES_IGNORE);

  spmvBufs_.rbuf_.resize(osetRecv);
  spmvBufs_.sbuf_.resize(osetSend);
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::spmv(const Cmpt* x, Cmpt* y) const {
  const int P = this->commcator_->getMySize();

  if (P != 1) {
    setupSpMVBuffers();
    for (ulabel i = 0; i < spmvBufs_.sind_.size(); i++)
      spmvBufs_.sbuf_[i] = x[spmvBufs_.sind_[i] - beginCol_];

    std::vector<MPI_Request> sreq(spmvBufs_.sranks_.size()),
        rreq(spmvBufs_.rranks_.size());

    exchange("S&R", spmvBufs_.sbuf_.data(), spmvBufs_.rbuf_.data());

#ifdef TIMERS
    TIMERS_START("SpMV diag part");
#endif

    Cmpt yrow;
    label nlocal, r;
    const label* ind;
    const Cmpt* val;
    const label* offd = offDiagStart_.data();
    const label* pt = this->ptr_;
    const label* ii = this->ind_;
    const Cmpt *v = this->val_, *xx;

#ifndef UNATPP

#ifdef OMP
#pragma omp parallel for
#endif
    for (r = 0; r < localRows_; r++) {
      nlocal = offd[r] - pt[r];
      ind = ii + pt[r];
      val = v + pt[r];
      xx = x - beginCol_;
      yrow = 0.;
      SparseDensePlusDot(yrow, xx, val, ind, nlocal);

      y[r] = yrow;
    }
#else
    if (localRows_ == localCols_) {
      if (topo_ == NULL) {
        const label diagNnz = localnnz_ - offDiagNnz_;
        topo_ = new unatpp::Topology(localRows_, localCols_, diagNnz);
        valDiag_ = new Cmpt[diagNnz];
        label count = 0;
        for (r = 0; r < localRows_; r++) {
          nlocal = offd[r] - pt[r];
          ind = ii + pt[r];
          val = v + pt[r];
          for (label j = 0; j < nlocal; j++) {
            topo_->add(r, ind[j] - beginCol_);
            valDiag_[count++] = val[j];
          }
        }
        ASSERT(count == diagNnz);

        topo_->sort_ents();
        topo_->switch_blksz(10);
        valTopo_ = new Cmpt[diagNnz];
        topo_->permute_matrix_data(valTopo_, valDiag_);
        delete[] valDiag_;
      }

      std::fill_n(y, localRows_, 0.0);
      iterateMV(
          *topo_,
          [] __host__ __device__(
              int e, int r, int c,
              unatpp::DataSet<unatpp::Matrix, unatpp::In, const Cmpt> val,
              unatpp::DataSet<unatpp::Row, unatpp::InOut, Cmpt> y,
              unatpp::DataSet<unatpp::Col, unatpp::In, const Cmpt> xx) {
            y[r] += xx[c] * val[e];
          },
          valTopo_, y, x);
    } else {
      for (r = 0; r < localRows_; r++) {
        nlocal = offd[r] - pt[r];
        ind = ii + pt[r];
        val = v + pt[r];
        xx = x - beginCol_;
        yrow = 0.;
        SparseDensePlusDot(yrow, xx, val, ind, nlocal);

        y[r] = yrow;
      }
    }

#endif
#ifdef TIMERS
    TIMERS_END("SpMV diag part");
#endif

    const label rowSize = spmvBufs_.rowIndex_.size();
    const label* compRow = spmvBufs_.compRow_.data();
    const label* prbuf = spmvBufs_.prbuf_.data();
    label* rowInd = spmvBufs_.rowIndex_.data();
    xx = spmvBufs_.rbuf_.data();

    this->commcator_->finishTask("S&R");

    for (label r = 0; r < rowSize; r++) {
      yrow = 0.;
      nlocal = compRow[r + 1] - compRow[r];
      val = v + offd[*rowInd];
      ind = prbuf + compRow[r];
      SparseDensePlusDot(yrow, xx, val, ind, nlocal);
      y[*rowInd++] += yrow;
    }

  } else {
    Cmpt yrow;
    label hij, r;
    const label* ind;
    const Cmpt* val;
    const label* pt = this->ptr_;
    const label* ii = this->ind_;
    const Cmpt* v = this->val_;
    const label m = this->gRows_;

#ifndef UNATPP
    for (r = 0; r < m; r++) {
      hij = pt[r + 1] - pt[r];
      ind = ii + pt[r];
      val = v + pt[r];
      yrow = 0.;
      SparseDensePlusDot(yrow, x, val, ind, hij);
      y[r] = yrow;
    }
#else
    if (localRows_ == localCols_) {
      if (topo_ == NULL) {
        const label diagNnz = localnnz_;
        topo_ = new unatpp::Topology(localRows_, localCols_, diagNnz);
        for (r = 0; r < localRows_; r++) {
          hij = pt[r + 1] - pt[r];
          ind = ii + pt[r];
          val = v + pt[r];
          for (label j = 0; j < hij; j++) {
            topo_->add(r, ind[j] - beginCol_);
          }
        }
        ASSERT(count == diagNnz);

        topo_->sort_ents();
        topo_->switch_blksz(10);
        valTopo_ = new Cmpt[diagNnz];
        topo_->permute_matrix_data(valTopo_, const_cast<Cmpt*>(v));
      }

      std::fill_n(y, localRows_, 0.0);
      iterateMV(
          *topo_,
          [] __host__ __device__(
              int e, int r, int c,
              unatpp::DataSet<unatpp::Matrix, unatpp::In, const Cmpt> val,
              unatpp::DataSet<unatpp::Row, unatpp::InOut, Cmpt> y,
              unatpp::DataSet<unatpp::Col, unatpp::In, const Cmpt> xx) {
            y[r] += xx[c] * val[e];
          },
          valTopo_, y, x);
    } else {
      for (r = 0; r < m; r++) {
        hij = pt[r + 1] - pt[r];
        ind = ii + pt[r];
        val = v + pt[r];
        yrow = 0.;
        SparseDensePlusDot(yrow, x, val, ind, hij);
        y[r] = yrow;
      }
    }
#endif
  }
}

template <typename Cmpt>
scalar CSRMatrix<Cmpt>::maxScaledRes(const Cmpt* x, const Cmpt* b) const {
  const int P = this->commcator_->getMySize();
  if (P != 1) {
    MPI_Comm comm = this->commcator_->getMpiComm();

    setupSpMVBuffers();

    for (ulabel i = 0; i < spmvBufs_.sind_.size(); i++)
      spmvBufs_.sbuf_[i] = x[spmvBufs_.sind_[i] - beginRow_];

    std::vector<MPI_Request> sreq(spmvBufs_.sranks_.size());
    for (ulabel p = 0; p < spmvBufs_.sranks_.size(); p++)
      MPI_Isend(spmvBufs_.sbuf_.data() + spmvBufs_.soff_[p],
                spmvBufs_.soff_[p + 1] - spmvBufs_.soff_[p], mpiType<Cmpt>(),
                spmvBufs_.sranks_[p], 0, comm, &sreq[p]);

    std::vector<MPI_Request> rreq(spmvBufs_.rranks_.size());
    for (ulabel p = 0; p < spmvBufs_.rranks_.size(); p++)
      MPI_Irecv(spmvBufs_.rbuf_.data() + spmvBufs_.roffs_[p],
                spmvBufs_.roffs_[p + 1] - spmvBufs_.roffs_[p], mpiType<Cmpt>(),
                spmvBufs_.rranks_[p], 0, comm, &rreq[p]);

    MPI_Waitall(rreq.size(), rreq.data(), MPI_STATUSES_IGNORE);

    scalar m = scalar(0.), tmp;
    std::vector<label>::iterator pbuf = spmvBufs_.prbuf_.begin();

    for (label r = 0; r < localRows_; r++) {
      Cmpt trueRes = b[r];
      scalar absRes = std::abs(b[r]);
      for (label j = this->ptr_[r]; j < offDiagStart_[r]; j++) {
        label c = this->ind_[j];
        trueRes -= this->val_[j] * x[c - beginRow_];
        absRes += std::abs(this->val_[j]) * std::abs(x[c - beginRow_]);
      }
      for (label j = offDiagStart_[r]; j < this->ptr_[r + 1]; j++) {
        trueRes -= this->val_[j] * spmvBufs_.rbuf_[*pbuf];
        absRes += std::abs(this->val_[j]) * std::abs(spmvBufs_.rbuf_[*pbuf]);
        pbuf++;
      }
      tmp = (absRes != 0) ? (std::abs(trueRes) / std::abs(absRes)) : 0.;
      m = std::max(m, tmp);
    }

    MPI_Waitall(sreq.size(), sreq.data(), MPI_STATUSES_IGNORE);
    MPI_Allreduce(MPI_IN_PLACE, &m, 1, mpiType<scalar>(), MPI_MAX, comm);
    return m;
  } else {
    scalar m = scalar(0.), tmp;
    for (label r = 0; r < this->gRows_; r++) {
      Cmpt trueRes = b[r];
      scalar absRes = std::abs(b[r]);
      const label hij = this->ptr_[r + 1];
      for (label j = this->ptr_[r]; j < hij; ++j) {
        const Cmpt v = this->val_[j];
        const label rj = this->ind_[j];
        trueRes -= v * x[rj];
        absRes += std::abs(v) * std::abs(x[rj]);
      }
      tmp = (absRes != 0) ? (std::abs(trueRes) / std::abs(absRes)) : 0.;
      m = std::max(m, tmp);
    }
    return m;
  }
}

template <typename Cmpt>
scalar CSRMatrix<Cmpt>::norm1() const {
  const int P = this->commcator_->getMySize();
  if (P != 1) {
    MPI_Comm comm = this->commcator_->getMpiComm();
    std::vector<scalar> n1(this->gRows_);
    for (label i = 0; i < localRows_; i++)
      for (label j = this->ptr_[i]; j < this->ptr_[i + 1]; j++)
        n1[this->ind_[j]] += std::abs(this->val_[j]);
    if (this->commcator_->getMyId() == 0)
      MPI_Reduce(MPI_IN_PLACE, n1.data(), n1.size(), mpiType<scalar>(), MPI_SUM,
                 0, comm);
    else
      MPI_Reduce(n1.data(), n1.data(), n1.size(), mpiType<scalar>(), MPI_SUM, 0,
                 comm);
    scalar nrm1 = 0;
    if (COMM::getGlobalId() == 0)
      nrm1 = *std::max_element(n1.begin(), n1.end());
    MPI_Bcast(&nrm1, 1, mpiType<scalar>(), 0, comm);
    return nrm1;
  } else {
    std::vector<scalar> n1(this->gRows_);
    for (label i = 0; i < this->gRows_; i++)
      for (label j = this->ptr_[i]; j < this->ptr_[i + 1]; j++)
        n1[this->ind_[j]] += std::abs(this->val_[j]);
    return *std::max_element(n1.begin(), n1.end());
  }
}

template <typename Cmpt>
CSRMatrix<Cmpt>* CSRMatrix<Cmpt>::gather() const {
  const int rank = this->commcator_->getMyId();
  const int P = this->commcator_->getMySize();
  MPI_Comm comm = this->commcator_->getMpiComm();
  if (rank == 0) {
    int* iwork = new int[2 * P];
    int* rcnts = iwork;
    int* displs = rcnts + P;
    for (label p = 0; p < P; p++) {
      rcnts[p] = vtxdist_[p + 1] - vtxdist_[p];
      displs[p] = vtxdist_[p];
    }
    CSRMatrix<Cmpt>* Aseq = new CSRMatrix<Cmpt>(this->gRows_, this->nnz_);
    MPI_Gatherv(this->ptr() + 1, localRows_, mpiType<label>(), Aseq->ptr() + 1,
                rcnts, displs, mpiType<label>(), 0, comm);
    Aseq->ptr(0) = 0;
    for (label p = 1; p < P; p++) {
      if (vtxdist_[p] > 0) {
        label p_start = Aseq->ptr(vtxdist_[p]);
        for (label r = vtxdist_[p]; r < vtxdist_[p + 1]; r++)
          Aseq->ptr(r + 1) += p_start;
      }
    }
    for (label p = 0; p < P; p++) {
      rcnts[p] = Aseq->ptr(vtxdist_[p + 1]) - Aseq->ptr(vtxdist_[p]);
      displs[p] = Aseq->ptr(vtxdist_[p]);
    }
    MPI_Gatherv(this->ind(), localnnz_, mpiType<label>(), Aseq->ind(), rcnts,
                displs, mpiType<label>(), 0, comm);
    MPI_Gatherv(this->val(), localnnz_, mpiType<Cmpt>(), Aseq->val(), rcnts,
                displs, mpiType<Cmpt>(), 0, comm);
    delete[] iwork;
    return Aseq;
  } else {
    MPI_Gatherv(this->ptr() + 1, localRows_, mpiType<label>(), NULL, NULL, NULL,
                mpiType<label>(), 0, comm);
    MPI_Gatherv(this->ind(), localnnz_, mpiType<label>(), NULL, NULL, NULL,
                mpiType<label>(), 0, comm);
    MPI_Gatherv(this->val(), localnnz_, mpiType<Cmpt>(), NULL, NULL, NULL,
                mpiType<Cmpt>(), 0, comm);
    return NULL;
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::getSymmPattern() {
  if (this->symm_) return;
  const int P = this->commcator_->getMySize();

  if (P == 1) {
    std::vector<label> a2Ctr(this->gRows_);

    for (label i = 0; i < this->gRows_; i++)
      a2Ctr[i] = this->ptr_[i + 1] - this->ptr_[i];

    bool change = false;

    for (label i = 0; i < this->gRows_; i++)
      for (label jj = this->ptr_[i]; jj < this->ptr_[i + 1]; jj++) {
        label kb = this->ptr_[this->ind_[jj]],
              ke = this->ptr_[this->ind_[jj] + 1];
        if (std::find(this->ind_ + kb, this->ind_ + ke, i) == this->ind_ + ke) {
          {
            a2Ctr[this->ind_[jj]]++;
            change = true;
          }
        }
      }
    if (change) {
      std::vector<label> a2Ptr(this->gRows_ + 1);
      a2Ptr[0] = 0;
      for (label i = 0; i < this->gRows_; i++)
        a2Ptr[i + 1] = a2Ptr[i] + a2Ctr[i];
      label newNnz = a2Ptr[this->gRows_] - a2Ptr[0];
      std::vector<label> a2Ind(newNnz);
      std::vector<Cmpt> a2Val(newNnz);
      this->nnz_ = newNnz;

      for (label i = 0; i < this->gRows_; i++) {
        a2Ctr[i] = a2Ptr[i] + this->ptr_[i + 1] - this->ptr_[i];
        for (label jj = this->ptr_[i], k = a2Ptr[i]; jj < this->ptr_[i + 1];
             jj++) {
          a2Ind[k] = this->ind_[jj];
          a2Val[k++] = this->val_[jj];
        }
      }

      for (label i = 0; i < this->gRows_; i++)
        for (label jj = this->ptr_[i]; jj < this->ptr_[i + 1]; jj++) {
          label kb = this->ptr_[this->ind_[jj]],
                ke = this->ptr_[this->ind_[jj] + 1];
          if (std::find(this->ind_ + kb, this->ind_ + ke, i) ==
              this->ind_ + ke) {
            label t = this->ind_[jj];

            {
              a2Ind[a2Ctr[t]] = i;
              a2Val[a2Ctr[t]] = Cmpt(0.);
              a2Ctr[t]++;
            }
          }
        }

      delete[] this->ind_;
      delete[] this->val_;
      this->ind_ = new label[newNnz];
      this->val_ = new Cmpt[newNnz];
      std::copy(a2Ptr.begin(), a2Ptr.end(), this->ptr_);
      std::copy(a2Ind.begin(), a2Ind.end(), this->ind_);
      std::copy(a2Val.begin(), a2Val.end(), this->val_);
    }
    this->symm_ = true;
  } else {
    MPI_Comm comm = this->commcator_->getMpiComm();
    std::vector<std::vector<IdxIJ> > sbuf(P);
    for (label r = 0; r < localRows_; r++)
      for (label j = offDiagStart_[r]; j < this->ptr_[r + 1]; j++) {
        label col = this->ind_[j];
        label row = r + beginRow_;
        label dest = std::upper_bound(vtxdist_.begin(), vtxdist_.end(), col) -
                     vtxdist_.begin() - 1;

        sbuf[dest].push_back(IdxIJ(col, row));
      }
    std::vector<IdxIJ> edges;
    allToallv(edges, comm, sbuf);
    std::stable_sort(edges.begin(), edges.end(), PATCH::compareIdxIJ);

    label* rowSums = new label[localRows_];
    for (label r = 0; r < localRows_; r++)
      rowSums[r] = this->ptr_[r + 1] - this->ptr_[r];
    label newNnz = localnnz_;
    std::vector<IdxIJ>::iterator ep = edges.begin();
    for (label r = 0; r < localRows_; r++) {
      while (ep != edges.end() && ep->i_ < r + beginRow_) ep++;
      if (ep == edges.end()) break;
      while (ep != edges.end() && ep->i_ == r + beginRow_) {
        label kb = offDiagStart_[r], ke = this->ptr_[r + 1];
        if (std::find(this->ind() + kb, this->ind() + ke, ep->j_) ==
            this->ind() + ke) {
          newNnz++;
          rowSums[r]++;
        }
        ep++;
      }
    }

    for (label r = 0; r < localRows_; r++)
      for (label j = this->ptr_[r]; j < offDiagStart_[r]; j++) {
        label lc = this->ind_[j] - beginRow_;
        label kb = this->ptr_[lc], ke = offDiagStart_[lc];
        if (std::find(this->ind() + kb, this->ind() + ke, r + beginRow_) ==
            this->ind() + ke) {
          rowSums[lc]++;
          newNnz++;
        }
      }
    if (newNnz != localnnz_) {
      localnnz_ = newNnz;

      std::vector<label> newPtr(localRows_ + 1);
      newPtr[0] = 0;
      for (label r = 0; r < localRows_; r++)
        newPtr[r + 1] = newPtr[r] + rowSums[r];
      std::vector<label> newInd(newNnz);
      std::vector<Cmpt> newVal(newNnz);

      for (label r = 0; r < localRows_; r++) {
        rowSums[r] = newPtr[r] + this->ptr_[r + 1] - this->ptr_[r];
        for (label j = this->ptr_[r], k = newPtr[r]; j < this->ptr_[r + 1];
             j++) {
          newInd[k] = this->ind_[j];
          newVal[k++] = this->val_[j];
        }
      }

      for (label r = 0; r < localRows_; r++)
        for (label j = this->ptr_[r]; j < offDiagStart_[r]; j++) {
          label lc = this->ind_[j] - beginRow_;
          label kb = this->ptr_[lc], ke = offDiagStart_[lc];
          if (std::find(this->ind() + kb, this->ind() + ke, r + beginRow_) ==
              this->ind() + ke) {
            newInd[rowSums[lc]] = r + beginRow_;
            newVal[rowSums[lc]] = Cmpt(0.);
            rowSums[lc]++;
          }
        }

      ep = edges.begin();
      for (label r = 0; r < localRows_; r++) {
        while (ep != edges.end() && ep->i_ < r + beginRow_) ep++;
        if (ep == edges.end()) break;
        while (ep != edges.end() && ep->i_ == r + beginRow_) {
          label kb = offDiagStart_[r], ke = this->ptr_[r + 1];
          if (std::find(this->ind() + kb, this->ind() + ke, ep->j_) ==
              this->ind() + ke) {
            newInd[rowSums[r]] = ep->j_;
            newVal[rowSums[r]] = Cmpt(0.);
            rowSums[r]++;
          }
          ep++;
        }
      }
      delete[] this->ind_;
      delete[] this->val_;
      this->ind_ = new label[newNnz];
      this->val_ = new Cmpt[newNnz];
      std::copy(newPtr.begin(), newPtr.end(), this->ptr_);
      std::copy(newInd.begin(), newInd.end(), this->ind_);
      std::copy(newVal.begin(), newVal.end(), this->val_);
    }
    delete[] rowSums;
    label totalNewNnz;
    this->commcator_->allReduce("sumLocalNewNnz", &newNnz, &totalNewNnz, 1,
                                COMM_LABEL, COMM_SUM);

    if (totalNewNnz != this->nnz_) {
      splitOffDiag();
      this->nnz_ = totalNewNnz;
    }
    this->symm_ = true;
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::spMV(Vector<Cmpt>& Apsi, const Vector<Cmpt>& psi) const {
  CSRMatrix<Cmpt>::spmv(psi.data(), Apsi.data());
}

template <typename Cmpt>
label CSRMatrix<Cmpt>::tsize() {
  label& gSize = SparseMatrix<Cmpt>::size();
  const label lSize = localRows_;
  if (gSize == 0) {
    this->commcator_->allReduce("sumLocalNewNnz", &lSize, &gSize, 1, COMM_LABEL,
                                COMM_SUM);
    this->commcator_->finishTask("sumLocalNewNnz");
  }
  return gSize;
}

template <typename Cmpt>
Vector<Cmpt>& CSRMatrix<Cmpt>::diag() const {
  if (diagPtr_ == NULL) {
    diagPtr_ = new Vector<Cmpt>(localRows_, this->commcator_);
    label nlocal, r, j;
    const label* ind;
    const Cmpt* val;
    const label* pt = this->ptr_;
    const label* ii = this->ind_;
    const Cmpt* v = this->val_;
    for (r = 0; r < localRows_; r++) {
      nlocal = pt[r + 1] - pt[r];
      ind = ii + pt[r];
      val = v + pt[r];
      for (j = 0; j < nlocal; j++) {
        if (r == ind[j] - beginRow_) (*diagPtr_)[r] = val[j];
      }
    }
  }
  return *diagPtr_;
}

template <typename Cmpt>
Interfaces<Cmpt>* CSRMatrix<Cmpt>::matrixInterfaces() const {
  COUT << "class CSRMatrix don't need class Interfaces" << ENDL;
  return NULL;
}
template <typename Cmpt>
void CSRMatrix<Cmpt>::initInterfaces(const Vector<Cmpt>& psi) const {
  COUT << "initInterfaces function for CSRMatrix did nothing" << ENDL;
}
template <typename Cmpt>
void CSRMatrix<Cmpt>::updateInterfaces(Vector<Cmpt>& Apsi) const {
  COUT << "updateInterfaces function for CSRMatrix did nothing" << ENDL;
}

template <typename Cmpt>
CSRMatrix<Cmpt>* CSRMatrix<Cmpt>::remoteRows(const CSRMatrix<Cmpt>& A,
                                             bool needValues) const {
  const SPMVBuffers<Cmpt>& spmvBufs = this->spmvBufs_;
  const label nrecv = spmvBufs.rranks_.size();
  const label nsend = spmvBufs.sranks_.size();
  std::vector<CSRMatrix<Cmpt> > sendRows(nsend);

  static const int tagPtr = 2001;
  static const int tagInd = 2002;
  static const int tagVal = 2003;
  std::vector<MPI_Request> sendPtrReq(nsend);
  std::vector<MPI_Request> sendIndReq(nsend);
  std::vector<MPI_Request> sendValReq(nsend);
  std::vector<MPI_Request> recvPtrReq(nrecv);
  std::vector<MPI_Request> recvIndReq(nrecv);
  std::vector<MPI_Request> recvValReq(nrecv);
  MPI_Comm mpiComm = this->commcator_->getMpiComm();

  for (label k = 0; k < nsend; k++) {
    CSRMatrix<Cmpt>& m = sendRows[k];
    const label beg = spmvBufs.soff_[k];
    const label end = spmvBufs.soff_[k + 1];
    m.lRows() = end - beg;
    m.ptr() = new label[m.lRows() + 1];
    m.lnnz() = 0;

    label ii = 0;

    for (label i = beg; i < end; ++i) {
      const label r = spmvBufs.sind_[i] - beginCol_;
      const label w = A.ptr(r + 1) - A.ptr(r);
      m.ptr(ii++) = w;
      m.lnnz() += w;
    }
    MPI_Isend(m.ptr(), m.lRows() * sizeof(label), MPI_BYTE,
              spmvBufs_.sranks_[k], tagPtr, mpiComm, &sendPtrReq[k]);

    m.ind() = new label[m.lnnz()];
    if (needValues) m.val() = new Cmpt[m.lnnz()];

    ii = 0;
    for (label i = beg; i < end; ++i) {
      const label r = spmvBufs.sind_[i] - beginCol_;
      for (label j = A.ptr(r); j < A.ptr(r + 1); j++) {
        m.ind(ii) = A.ind(j);
        if (needValues) m.val(ii) = A.val(j);
        ++ii;
      }
    }
    MPI_Isend(m.ind(), (m.lnnz()) * sizeof(label), MPI_BYTE,
              spmvBufs_.sranks_[k], tagInd, mpiComm, &sendIndReq[k]);

    if (needValues)
      MPI_Isend(m.val(), (m.lnnz()) * sizeof(Cmpt), MPI_BYTE,
                spmvBufs_.sranks_[k], tagVal, mpiComm, &sendValReq[k]);
  }

  CSRMatrix<Cmpt>* S = new CSRMatrix<Cmpt>(A.comm());
  S->lRows() = spmvBufs.rind_.size();
  S->ptr() = new label[S->lRows() + 1];
  S->ptr(0) = 0;

  for (label k = 0; k < nrecv; k++) {
    const label beg = spmvBufs.roffs_[k];
    const label end = spmvBufs.roffs_[k + 1];
    MPI_Irecv(&S->ptr(beg + 1), (end - beg) * sizeof(label), MPI_BYTE,
              spmvBufs_.rranks_[k], tagPtr, mpiComm, &recvPtrReq[k]);
  }
  MPI_Waitall(recvPtrReq.size(), recvPtrReq.data(), MPI_STATUSES_IGNORE);

  std::partial_sum(S->ptr(), S->ptr() + S->lRows() + 1, S->ptr());
  S->lnnz() = S->ptr(S->lRows());
  S->ind() = new label[S->lnnz()];
  if (needValues) S->val() = new Cmpt[S->lnnz()];
  for (label k = 0; k < nrecv; k++) {
    const label beg = spmvBufs.roffs_[k];
    const label end = spmvBufs.roffs_[k + 1];
    const label cbeg = S->ptr(beg);
    const label cend = S->ptr(end);
    MPI_Irecv(&S->ind(cbeg), (cend - cbeg) * sizeof(label), MPI_BYTE,
              spmvBufs_.rranks_[k], tagInd, mpiComm, &recvIndReq[k]);

    if (needValues)
      MPI_Irecv(&S->val(cbeg), (cend - cbeg) * sizeof(Cmpt), MPI_BYTE,
                spmvBufs_.rranks_[k], tagVal, mpiComm, &recvValReq[k]);
  }

  MPI_Waitall(sendPtrReq.size(), sendPtrReq.data(), MPI_STATUSES_IGNORE);
  MPI_Waitall(sendIndReq.size(), sendIndReq.data(), MPI_STATUSES_IGNORE);
  MPI_Waitall(recvIndReq.size(), recvIndReq.data(), MPI_STATUSES_IGNORE);

  if (needValues) {
    MPI_Waitall(sendValReq.size(), sendValReq.data(), MPI_STATUSES_IGNORE);
    MPI_Waitall(recvValReq.size(), recvValReq.data(), MPI_STATUSES_IGNORE);
  }

  return S;
}

template <typename Cmpt>
csr<Cmpt>* CSRMatrix<Cmpt>::extractDiagPart() const {
  csr<Cmpt>* t = new csr<Cmpt>();
  csr<Cmpt>& T = *t;
  T.setSize(localRows_, localCols_);
  const label nnz = localnnz_ - offDiagNnz_;
  for (label i = 0; i < localRows_; ++i) {
    T.ptr[i + 1] = offDiagStart_[i] - this->ptr_[i];
  }
  T.setNnz(T.scanRowSize());
  ASSERT(T.nnz == nnz);
  label count = 0;
  for (label i = 0; i < localRows_; ++i) {
    for (label j = this->ptr_[i], e = offDiagStart_[i]; j < e; ++j) {
      T.col[count] = this->ind_[j] - beginCol_;
      T.val[count++] = this->val_[j];
    }
  }

  return t;
}

template <typename Cmpt>
csr<Cmpt>* CSRMatrix<Cmpt>::extractOffdiagPart() const {
  csr<Cmpt>* t = new csr<Cmpt>();
  csr<Cmpt>& T = *t;
  unordered_map<label, std::pair<label, label> >& g2l = spmvBufs_.rind_;
  T.setSize(localRows_, spmvBufs_.rbuf_.size());
  for (label i = 0; i < localRows_; ++i) {
    T.ptr[i + 1] = this->ptr_[i + 1] - offDiagStart_[i];
  }
  T.setNnz(T.scanRowSize());
  ASSERT(T.nnz == offDiagNnz_);
  label count = 0;
  for (label i = 0; i < localRows_; ++i) {
    for (label j = offDiagStart_[i], e = this->ptr_[i + 1]; j < e; ++j) {
      T.col[count] = g2l[this->ind_[j]].second;
      T.val[count++] = this->val_[j];
    }
  }

  return t;
}

template <typename Cmpt>
CSRMatrix<Cmpt>* CSRMatrix<Cmpt>::transpose() const {
  Communicator* comm = this->commcator_;
  MPI_Comm mpiComm = comm->getMpiComm();
  const label nProcs = comm->getMySize();
  if (nProcs > 1) {
    SPMVBuffers<Cmpt>& spmvbuf = this->spMVBuffers();
    const label rsize = spmvbuf.sranks_.size();
    const label ssize = spmvbuf.rranks_.size();
    static const int tagCnt = 2001;
    static const int tagCol = 2002;
    static const int tagVal = 2003;
    std::vector<MPI_Request> recvCntReq(rsize);
    std::vector<MPI_Request> recvColReq(rsize);
    std::vector<MPI_Request> recvValReq(rsize);
    std::vector<MPI_Request> sendCntReq(ssize);
    std::vector<MPI_Request> sendColReq(ssize);
    std::vector<MPI_Request> sendValReq(ssize);
    csr<Cmpt>* Aloc = extractDiagPart();
    csr<Cmpt>* Arem = extractOffdiagPart();
    csr<Cmpt>* trem = Arem->transpose();
    csr<Cmpt>& Trem = *trem;
    delete Arem;

    label locBeg = vtxdist_[comm->getMyId()];
    for (label i = 0; i < Trem.nnz; ++i) Trem.col[i] += locBeg;

    std::vector<label> rowSize(Trem.nrows);
    for (label i = 0; i < Trem.nrows; ++i)
      rowSize[i] = Trem.ptr[i + 1] - Trem.ptr[i];

    std::vector<label> remPtr(spmvBufs_.sind_.size() + 1);
    remPtr[0] = 0;

    for (label i = 0; i < rsize; ++i) {
      label beg = spmvBufs_.soff_[i];
      label end = spmvBufs_.soff_[i + 1];

      MPI_Irecv(&remPtr[beg + 1], (end - beg) * sizeof(label), MPI_BYTE,
                spmvBufs_.sranks_[i], tagCnt, mpiComm, &recvCntReq[i]);
    }

    for (label i = 0; i < ssize; ++i) {
      label beg = spmvBufs_.roffs_[i];
      label end = spmvBufs_.roffs_[i + 1];

      MPI_Isend(&rowSize[beg], (end - beg) * sizeof(label), MPI_BYTE,
                spmvBufs_.rranks_[i], tagCnt, mpiComm, &sendCntReq[i]);
    }

    MPI_Waitall(recvCntReq.size(), recvCntReq.data(), MPI_STATUSES_IGNORE);

    std::partial_sum(remPtr.begin(), remPtr.end(), remPtr.begin());

    std::vector<label> remCol(remPtr.back());
    std::vector<Cmpt> remVal(remPtr.back());

    for (label i = 0; i < rsize; ++i) {
      label rbeg = spmvBufs_.soff_[i];
      label rend = spmvBufs_.soff_[i + 1];

      label cbeg = remPtr[rbeg];
      label cend = remPtr[rend];

      MPI_Irecv(&remCol[cbeg], (cend - cbeg) * sizeof(label), MPI_BYTE,
                spmvBufs_.sranks_[i], tagCol, mpiComm, &recvColReq[i]);

      MPI_Irecv(&remVal[cbeg], (cend - cbeg) * sizeof(Cmpt), MPI_BYTE,
                spmvBufs_.sranks_[i], tagVal, mpiComm, &recvValReq[i]);
    }

    for (label i = 0; i < ssize; ++i) {
      label rbeg = spmvBufs_.roffs_[i];
      label rend = spmvBufs_.roffs_[i + 1];

      label cbeg = Trem.ptr[rbeg];
      label cend = Trem.ptr[rend];

      MPI_Isend(&Trem.col[cbeg], (cend - cbeg) * sizeof(label), MPI_BYTE,
                spmvBufs_.rranks_[i], tagCol, mpiComm, &sendColReq[i]);

      MPI_Isend(&Trem.val[cbeg], (cend - cbeg) * sizeof(Cmpt), MPI_BYTE,
                spmvBufs_.rranks_[i], tagVal, mpiComm, &sendValReq[i]);
    }

    csr<Cmpt>* tremt = new csr<Cmpt>();
    csr<Cmpt>& Tremt = *tremt;
    Tremt.setSize(localCols_, 0);

    for (ulabel i = 0; i < spmvBufs_.sind_.size(); ++i)
      Tremt.ptr[1 + spmvBufs_.sind_[i] - beginCol_] +=
          remPtr[i + 1] - remPtr[i];

    Tremt.setNnz(Tremt.scanRowSize());

    csr<Cmpt>* tloc = Aloc->transpose();
    csr<Cmpt>& Tloc = *tloc;
    delete Aloc;

    MPI_Waitall(recvColReq.size(), recvColReq.data(), MPI_STATUSES_IGNORE);
    MPI_Waitall(recvValReq.size(), recvValReq.data(), MPI_STATUSES_IGNORE);

    for (ulabel i = 0; i < spmvBufs_.sind_.size(); ++i) {
      label row = spmvBufs_.sind_[i] - beginCol_;
      label head = Tremt.ptr[row];

      for (label j = remPtr[i]; j < remPtr[i + 1]; ++j, ++head) {
        Tremt.col[head] = remCol[j];
        Tremt.val[head] = remVal[j];
      }

      Tremt.ptr[row] = head;
    }

    std::rotate(Tremt.ptr, Tremt.ptr + localCols_, Tremt.ptr + localCols_ + 1);
    Tremt.ptr[0] = 0;

    MPI_Waitall(sendCntReq.size(), sendCntReq.data(), MPI_STATUSES_IGNORE);
    MPI_Waitall(sendColReq.size(), sendColReq.data(), MPI_STATUSES_IGNORE);
    MPI_Waitall(sendValReq.size(), sendValReq.data(), MPI_STATUSES_IGNORE);
    delete trem;

    CSRMatrix<Cmpt>* t = new CSRMatrix<Cmpt>(comm);
    CSRMatrix<Cmpt>& T = *t;
    T.lRows() = localCols_;
    T.offDiagStart().resize(localCols_);
    T.lCols() = localRows_;
    T.beginRow() = beginCol_;
    T.beginCol() = beginRow_;
    T.lnnz() = Tloc.nnz + Tremt.nnz;
    T.offDiagNnz() = Tremt.nnz;
    T.vtxdist() = vtxdistCol_;
    T.vtxdistCol() = vtxdist_;
    T.ptr() = new label[localCols_ + 1];
    T.ind() = new label[T.lnnz()];
    T.val() = new Cmpt[T.lnnz()];
    T.setSymm(this->symm());
    label* const Tptr = T.ptr();
    label* const Tind = T.ind();
    Cmpt* const Tval = T.val();
    label* const Toffdiag = T.offDiagStart().data();

    Tptr[0] = 0;
    for (label i = 0; i < localCols_; ++i) {
      Tptr[i] = Tloc.ptr[i] + Tremt.ptr[i];
      Toffdiag[i] = Tloc.ptr[i + 1] + Tremt.ptr[i];
    }
    Tptr[localCols_] = Tloc.ptr[localCols_] + Tremt.ptr[localCols_];

    label count = 0;
    for (label i = 0; i < localCols_; ++i) {
      for (label j = Tloc.ptr[i]; j < Tloc.ptr[i + 1]; ++j) {
        const label inda = Tloc.col[j];
        const Cmpt ja = Tloc.val[j];
        Tind[count] = inda + beginRow_;
        Tval[count++] = ja;
      }
      for (label j = Tremt.ptr[i]; j < Tremt.ptr[i + 1]; ++j) {
        const label inda = Tremt.col[j];
        const Cmpt ja = Tremt.val[j];
        Tind[count] = inda;
        Tval[count++] = ja;
      }
    }

    delete tloc;
    delete tremt;
    return t;
  } else {
    csr<Cmpt>* Aloc = extractDiagPart();
    csr<Cmpt>* tloc = Aloc->transpose();
    csr<Cmpt>& Tloc = *tloc;
    delete Aloc;

    CSRMatrix<Cmpt>* t = new CSRMatrix<Cmpt>(comm);
    CSRMatrix<Cmpt>& T = *t;
    T.lRows() = localCols_;
    T.offDiagStart().resize(localCols_);
    T.lCols() = localRows_;
    T.beginRow() = beginCol_;
    T.beginCol() = beginRow_;
    T.rows() = this->cols();
    T.cols() = this->rows();
    T.lnnz() = Tloc.nnz;
    T.offDiagNnz() = 0;
    T.vtxdist() = vtxdistCol_;
    T.vtxdistCol() = vtxdist_;
    T.ptr() = new label[localCols_ + 1];
    T.ind() = new label[T.lnnz()];
    T.val() = new Cmpt[T.lnnz()];
    T.setSymm(this->symm());
    label* const Tptr = T.ptr();
    label* const Tind = T.ind();
    Cmpt* const Tval = T.val();

    Tptr[0] = 0;
    for (label i = 0; i < localCols_; ++i) {
      Tptr[i] = Tloc.ptr[i];
    }
    Tptr[localCols_] = Tloc.ptr[localCols_];

    label count = 0;
    for (label i = 0; i < localCols_; ++i) {
      for (label j = Tloc.ptr[i]; j < Tloc.ptr[i + 1]; ++j) {
        const label inda = Tloc.col[j];
        const Cmpt ja = Tloc.val[j];
        Tind[count] = inda + beginRow_;
        Tval[count++] = ja;
      }
    }

    delete tloc;
    return t;
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::scale(scalar s) {
  if (s == 1.) return;
  const label n = localRows_;
  const label* ptr = this->ptr_;
  Cmpt* val = this->val_;
  for (label i = 0; i < n; ++i) {
    for (label j = ptr[i]; j < ptr[i + 1]; ++j) {
      val[j] *= s;
    }
  }
}

template <typename Cmpt>
static Array<Array<Cmpt> >* val2D_;

template <typename Cmpt>
void CSRMatrix<Cmpt>::setMatrixTopology(const label rowSize,
                                        const labelVector& upperAddr,
                                        const labelVector& lowerAddr,
                                        const bool reUse) {
  const int nProcs = this->commcator_->getMySize();
  vtxdist_.resize(nProcs + 1);
  localRows_ = rowSize;
  localCols_ = localRows_;
  this->commcator_->allGather("getGlobalDimCSR", &localRows_, sizeof(label),
                              vtxdist_.data() + 1, sizeof(label));
  vtxdist_[0] = 0;
  lduConstrctor_.upper_ = new labelVector(upperAddr);
  lduConstrctor_.lower_ = new labelVector(lowerAddr);

  localnnz_ = rowSize + 2 * upperAddr.size();
  this->nnz_ = -1;
  this->symm_ = (&upperAddr == &lowerAddr);
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::createInterfacesTopology(const int nNeiProcs,
                                               const int* destRank,
                                               const label* offDiagRows,
                                               const label* offDiagCols,
                                               const label* offDiagStarts) {
  const int myId = this->commcator_->getMyId();
  const int nProcs = this->commcator_->getMySize();
  this->commcator_->finishTask("getGlobalDimCSR");
  for (int i = 0; i < nProcs; i++) vtxdist_[i + 1] += vtxdist_[i];
  beginRow_ = vtxdist_[myId];
  beginCol_ = beginRow_;
  vtxdistCol_ = vtxdist_;
  this->gRows_ = vtxdist_[nProcs];
  this->gCols_ = this->gRows_;

  offDiagNnz_ = 0;
  forAll(i, nNeiProcs) offDiagNnz_ += offDiagStarts[i + 1] - offDiagStarts[i];
  localnnz_ += offDiagNnz_;

  this->ptr_ = new label[localRows_ + 1];
  this->ind_ = new label[localnnz_];
  this->val_ = new Cmpt[localnnz_];

  Array<Array<label> > ind2D(localRows_);
  offDiagStart_.resize(localRows_);

  for (label i = 0; i < localRows_; i++) {
    ind2D[i].reserve(localnnz_ / localRows_);
  }

  const label faceNum = lduConstrctor_.upper_->size();
  const label* upperAddr = lduConstrctor_.upper_->data();
  const label* lowerAddr = lduConstrctor_.lower_->data();
  for (label i = 0; i < faceNum; i++) {
    ind2D[upperAddr[i]].push_back(lowerAddr[i] + beginRow_);
  }
  for (label i = 0; i < localRows_; i++) {
    ind2D[i].push_back(i + beginRow_);
  }
  for (label i = 0; i < faceNum; i++) {
    ind2D[lowerAddr[i]].push_back(upperAddr[i] + beginRow_);
  }
  for (label i = 0; i < localRows_; i++) offDiagStart_[i] = ind2D[i].size();

  Array<label> idx(offDiagNnz_);
  UNAP::iota(idx.begin(), idx.end(), 0);
  Comp compare(offDiagRows, offDiagCols);
  std::stable_sort(idx.begin(), idx.end(), compare);
  lduConstrctor_.order_ = new labelVector(idx);
  forAll(i, offDiagNnz_) {
    ind2D[offDiagRows[idx[i]]].push_back(offDiagCols[idx[i]]);
  }

  lduConstrctor_.offDiagRows_ = new Vector<label>(offDiagRows, offDiagNnz_);

  this->ptr_[0] = 0;
  forAll(i, localRows_) {
    this->ptr_[i + 1] = this->ptr_[i] + ind2D[i].size();
    offDiagStart_[i] += this->ptr_[i];

    std::copy(ind2D[i].begin(), ind2D[i].end(), this->ind_ + this->ptr_[i]);
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::setMatrixCoeffients(const Vector<Cmpt>& diag,
                                          const Vector<Cmpt>& upper,
                                          const Vector<Cmpt>& lower,
                                          const bool reuse) {
  DELETE_POINTER(diagPtr_);
  diagPtr_ = new Vector<Cmpt>(diag);
  val2D_<Cmpt> = new Array<Array<Cmpt> >(localRows_);

  Array<Array<Cmpt> >& val2D = *val2D_<Cmpt>;

  for (label i = 0; i < localRows_; i++) {
    val2D[i].reserve(localnnz_ / localRows_);
  }

  const label faceNum = lduConstrctor_.upper_->size();
  const label* upperAddr = lduConstrctor_.upper_->data();
  const label* lowerAddr = lduConstrctor_.lower_->data();
  const Cmpt* upperP = upper.data();
  const Cmpt* lowerP = lower.data();
  for (label i = 0; i < faceNum; i++) {
    val2D[upperAddr[i]].push_back(lowerP[i]);
  }
  for (label i = 0; i < localRows_; i++) {
    val2D[i].push_back((*diagPtr_)[i]);
  }
  for (label i = 0; i < faceNum; i++) {
    val2D[lowerAddr[i]].push_back(upperP[i]);
  }
}

template <typename Cmpt>
void CSRMatrix<Cmpt>::fillInterfacesCofficients(const label* offDiagStarts,
                                                const Cmpt* offDiagCoeffs) {
  Array<Array<Cmpt> >& val2D = *val2D_<Cmpt>;
  const label* idx = lduConstrctor_.order_->data();
  const label* rtmp = lduConstrctor_.offDiagRows_->data();
  forAll(i, offDiagNnz_) {
    val2D[rtmp[idx[i]]].push_back(offDiagCoeffs[idx[i]]);
  }
  forAll(i, localRows_) {
    std::copy(val2D[i].begin(), val2D[i].end(), this->val_ + this->ptr_[i]);
  }
  DELETE_POINTER(val2D_<Cmpt>)
}

template <>
scalar CSRMatrix<scalar>::spectralR(label iters) const {
  scalar radius = 0;
  if (iters <= 0) {
#ifdef OMP
#pragma omp parallel
#endif
    {
      scalar emax = 0;
      scalar dia = 1.0;
#ifdef OMP
#pragma omp for nowait
#endif
      for (label i = 0; i < lRows(); ++i) {
        scalar s = 0;
        for (label j = this->ptr_[i], e = this->ptr_[i + 1]; j < e; ++j) {
          label c = this->ind_[j];
          scalar v = this->val_[j];

          s += std::abs(v);

          if (c == i + beginRow_) dia = v;
        }

        s *= std::abs(1. / dia);
        emax = std::max(emax, s);
      }

#ifdef OMP
#pragma omp critical
#endif
      radius = std::max(radius, emax);
    }
  } else {
  }

  return radius < 0 ? static_cast<scalar>(2) : radius;
}

template <typename Cmpt>
CSRMatrix<Cmpt>* product(const CSRMatrix<Cmpt>& A, const CSRMatrix<Cmpt>& B) {
  const label Alrows = A.lRows();
  const label Alcols = A.lCols();
  const label Abegin = A.beginCol();
  const label Blrows = B.lRows();
  const label Blcols = B.lCols();
  const label Bbegin = B.beginCol();
  const label Bend = Bbegin + Blcols;
  ASSERT(Alcols == Blrows);
  const label *Aptr = A.ptr(), *Aind = A.ind(), *Bptr = B.ptr(),
              *Bind = B.ind();
  const int nProcs = A.comm()->getMySize();
  if (nProcs > 1) {
    SPMVBuffers<Cmpt>& Aspmvbuf = A.spMVBuffers();
    B.spMVBuffers();

    const label* Aoffdiag = A.offDiagStart().data();
    const label* Boffdiag = B.offDiagStart().data();
    const Cmpt *Aval = A.val(), *Bval = B.val();

    CSRMatrix<Cmpt>* b_nbr = A.remoteRows(B);
    CSRMatrix<Cmpt>& B_nbr = *b_nbr;

    CSRMatrix<Cmpt>* C = new CSRMatrix<Cmpt>(A.comm());
    C->lRows() = Alrows, C->lCols() = Blcols;
    C->beginRow() = A.beginRow(), C->beginCol() = Bbegin;
    C->ptr() = new label[Alrows + 1];
    C->vtxdist() = A.vtxdist();
    C->vtxdistCol() = B.vtxdistCol();
    std::vector<label>& Coffdiag = C->offDiagStart();
    Coffdiag.resize(Alrows);
    label* const Cptr = C->ptr();
    Cptr[0] = 0;

    std::vector<label> remCols;
    label ind = B_nbr.lnnz();
    remCols.resize(ind + B.offDiagNnz());
    std::copy(B_nbr.ind(), B_nbr.ind() + ind, remCols.begin());
    for (label i = 0; i < Blrows; i++) {
      for (label j = Boffdiag[i]; j < Bptr[i + 1]; j++) {
        remCols[ind++] = Bind[j];
      }
    }
    std::sort(remCols.begin(), remCols.end());
    remCols.erase(std::unique(remCols.begin(), remCols.end()), remCols.end());
    label count = 0;
    unordered_map<label, label> remIdx;
    for (ulabel i = 0; i < remCols.size(); i++) {
      const label c = remCols[i];
      if (c >= Bbegin && c < Bend) continue;
      remIdx[c] = count++;
    }

#ifdef OMP
#pragma omp parallel
#endif

    {
      std::vector<label> locMarker(Blcols, -1);
      std::vector<label> remMarker(count, -1);
      label offDiagNnzCount = 0;
#ifdef OMP
#pragma omp for
#endif
      for (label ia = 0; ia < Alrows; ++ia) {
        label locCols = 0, remCols = 0;

        for (label ja = Aptr[ia]; ja < Aoffdiag[ia]; ja++) {
          const label ca = Aind[ja] - Abegin;
          for (label jb = Bptr[ca]; jb < Bptr[ca + 1]; jb++) {
            label cb = Bind[jb];
            if (cb >= Bbegin && cb < Bend) {
              cb -= Bbegin;
              if (locMarker[cb] != ia) {
                locMarker[cb] = ia;
                ++locCols;
              }
            } else {
              cb = remIdx[cb];
              if (remMarker[cb] != ia) {
                remMarker[cb] = ia;
                ++remCols;
              }
            }
          }
        }

        for (label ja = Aoffdiag[ia]; ja < Aptr[ia + 1]; ja++) {
          const label ca = Aspmvbuf.rind_[Aind[ja]].second;
          for (label jb = B_nbr.ptr(ca); jb < B_nbr.ptr(ca + 1); jb++) {
            label cb = B_nbr.ind(jb);
            if (cb >= Bbegin && cb < Bend) {
              cb -= Bbegin;
              if (locMarker[cb] != ia) {
                locMarker[cb] = ia;
                ++locCols;
              }
            } else {
              cb = remIdx[cb];
              if (remMarker[cb] != ia) {
                remMarker[cb] = ia;
                ++remCols;
              }
            }
          }
        }
        Coffdiag[ia] = locCols;
        Cptr[ia + 1] = locCols + remCols;
        offDiagNnzCount += remCols;
      }
      C->offDiagNnz() = offDiagNnzCount;
    }

    std::partial_sum(Cptr, Cptr + Alrows + 1, Cptr);
    for (label ia = 0; ia < Alrows; ia++) Coffdiag[ia] += Cptr[ia];
    C->lnnz() = Cptr[Alrows];
    C->ind() = new label[C->lnnz()];
    C->val() = new Cmpt[C->lnnz()];
    label* Cind = C->ind();
    Cmpt* Cval = C->val();

#ifdef OMP
#pragma omp parallel
#endif

    {
      std::vector<label> locMarker(Blcols, -1);
      std::vector<label> remMarker(count, -1);
#ifdef OMP
#pragma omp for
#endif
      for (label ia = 0; ia < Alrows; ++ia) {
        const label locBegin = Cptr[ia];
        const label remBegin = Coffdiag[ia];
        label locEnd = locBegin;
        label remEnd = remBegin;

        for (label ja = Aptr[ia]; ja < Aoffdiag[ia]; ja++) {
          const label ca = Aind[ja] - Abegin;
          const Cmpt va = Aval[ja];
          for (label jb = Bptr[ca]; jb < Bptr[ca + 1]; jb++) {
            label cb = Bind[jb];
            const Cmpt vb = Bval[jb];
            if (cb >= Bbegin && cb < Bend) {
              cb -= Bbegin;
              if (locMarker[cb] < locBegin) {
                locMarker[cb] = locEnd;

                Cind[locEnd] = cb + Bbegin;
                Cval[locEnd] = va * vb;
                ++locEnd;
              } else {
                Cval[locMarker[cb]] += va * vb;
              }
            } else {
              const label gb = remIdx[cb];
              if (remMarker[gb] < remBegin) {
                remMarker[gb] = remEnd;

                Cind[remEnd] = cb;
                Cval[remEnd] = va * vb;
                ++remEnd;
              } else {
                Cval[remMarker[gb]] += va * vb;
              }
            }
          }
        }

        for (label ja = Aoffdiag[ia]; ja < Aptr[ia + 1]; ja++) {
          const label ca = Aspmvbuf.rind_[Aind[ja]].second;
          const Cmpt va = Aval[ja];
          for (label jb = B_nbr.ptr(ca); jb < B_nbr.ptr(ca + 1); jb++) {
            label cb = B_nbr.ind(jb);
            const Cmpt vb = B_nbr.val(jb);
            if (cb >= Bbegin && cb < Bend) {
              cb -= Bbegin;
              if (locMarker[cb] < locBegin) {
                locMarker[cb] = locEnd;

                Cind[locEnd] = cb + Bbegin;
                Cval[locEnd] = va * vb;
                ++locEnd;
              } else {
                Cval[locMarker[cb]] += va * vb;
              }
            } else {
              const label gb = remIdx[cb];
              if (remMarker[gb] < remBegin) {
                remMarker[gb] = remEnd;

                Cind[remEnd] = cb;
                Cval[remEnd] = va * vb;
                ++remEnd;
              } else {
                Cval[remMarker[gb]] += va * vb;
              }
            }
          }
        }
      }
    }
    delete b_nbr;
    return C;
  } else {
    const Cmpt *Aval = A.val(), *Bval = B.val();

    CSRMatrix<Cmpt>* C = new CSRMatrix<Cmpt>(A.comm());
    C->lRows() = Alrows, C->lCols() = Blcols;
    C->beginRow() = A.beginRow(), C->beginCol() = Bbegin;
    C->rows() = Alrows, C->cols() = Blcols;
    C->ptr() = new label[Alrows + 1];
    C->vtxdist() = A.vtxdist();
    C->vtxdistCol() = B.vtxdistCol();
    std::vector<label>& Coffdiag = C->offDiagStart();
    Coffdiag.resize(Alrows);
    label* const Cptr = C->ptr();
    Cptr[0] = 0;

#ifdef OMP
#pragma omp parallel
#endif

    {
      std::vector<label> locMarker(Blcols, -1);
#ifdef OMP
#pragma omp for
#endif
      for (label ia = 0; ia < Alrows; ++ia) {
        label locCols = 0;

        for (label ja = Aptr[ia]; ja < Aptr[ia + 1]; ja++) {
          const label ca = Aind[ja] - Abegin;
          for (label jb = Bptr[ca]; jb < Bptr[ca + 1]; jb++) {
            label cb = Bind[jb];
            cb -= Bbegin;
            if (locMarker[cb] != ia) {
              locMarker[cb] = ia;
              ++locCols;
            }
          }
        }

        Coffdiag[ia] = locCols;
        Cptr[ia + 1] = locCols;
      }
    }

    std::partial_sum(Cptr, Cptr + Alrows + 1, Cptr);
    for (label ia = 0; ia < Alrows; ia++) Coffdiag[ia] += Cptr[ia];
    C->lnnz() = Cptr[Alrows];
    C->ind() = new label[C->lnnz()];
    C->val() = new Cmpt[C->lnnz()];
    label* Cind = C->ind();
    Cmpt* Cval = C->val();

#ifdef OMP
#pragma omp parallel
#endif

    {
      std::vector<label> locMarker(Blcols, -1);
#ifdef OMP
#pragma omp for
#endif
      for (label ia = 0; ia < Alrows; ++ia) {
        const label locBegin = Cptr[ia];
        label locEnd = locBegin;

        for (label ja = Aptr[ia]; ja < Aptr[ia + 1]; ja++) {
          const label ca = Aind[ja] - Abegin;
          const Cmpt va = Aval[ja];
          for (label jb = Bptr[ca]; jb < Bptr[ca + 1]; jb++) {
            label cb = Bind[jb];
            const Cmpt vb = Bval[jb];
            cb -= Bbegin;
            if (locMarker[cb] < locBegin) {
              locMarker[cb] = locEnd;

              Cind[locEnd] = cb + Bbegin;
              Cval[locEnd] = va * vb;
              ++locEnd;
            } else {
              Cval[locMarker[cb]] += va * vb;
            }
          }
        }
      }
    }
    return C;
  }
}

template <typename Cmpt>
void readMatrixMarket(CSRMatrix<Cmpt>& A, const char* filename,
                      bool zeroBased) {
  enum MMsym { GENERAL, SYMMETRIC, SKEWSYMMETRIC, HERMITIAN };
  FILE* fp = fopen(filename, "r");
  if (fp == NULL) {
    std::cerr << "ERROR: could not read file";
    exit(1);
  }
  const label maxCline = 256;
  char cline[maxCline];
  if (fgets(cline, maxCline, fp) == NULL) {
    std::cerr << "ERROR: could not read from file" << std::endl;
    exit(1);
  }
  if (strstr(cline, "pattern")) {
    std::cerr << "ERROR: This is not a matrix,"
              << " but just a sparsity pattern" << std::endl;
    exit(1);
  } else if (strstr(cline, "complex")) {
    if (!isComplex<Cmpt>()) {
      fclose(fp);
      throw "ERROR: Complex matrix";
    }
  }
  MMsym s = GENERAL;
  if (strstr(cline, "skew-symmetric")) {
    s = SKEWSYMMETRIC;
    A.setSymm(true);
  } else if (strstr(cline, "symmetric")) {
    s = SYMMETRIC;
    A.setSymm(true);
  } else if (strstr(cline, "hermitian")) {
    s = HERMITIAN;
    A.setSymm(true);
  } else {
    A.setSymm(false);
  }

  while (fgets(cline, maxCline, fp)) {
    if (cline[0] != '%') {
      label m, in, innz;
#ifdef LABEL_INT64
      sscanf(cline, "%ld %ld %ld", &m, &in, &innz);
#else
      sscanf(cline, "%d %d %d", &m, &in, &innz);
#endif
      A.nnz() = static_cast<label>(innz);
      A.SparseMatrix<Cmpt>::size() = static_cast<label>(in);
      if (s != GENERAL) A.nnz() *= 2;
      if (m != A.tsize()) {
        std::cerr << "ERROR: matrix is not square!" << std::endl;
        exit(1);
      }
      break;
    }
  }

  const int nProcs = A.comm()->getMySize();
  const int myId = A.comm()->getMyId();

  label avgColNum = A.tsize() / nProcs;
  label remainder = A.tsize() % nProcs;
  label myBeginInd, myEndInd;

  if (myId < remainder) {
    myBeginInd = myId * (avgColNum + 1);
    myEndInd = (myId + 1) * (avgColNum + 1);

  } else {
    myBeginInd = myId * avgColNum + remainder;
    myEndInd = (myId + 1) * avgColNum + remainder;
  }

  std::vector<Triplet<Cmpt> > Atuple;
  Atuple.reserve(A.nnz() / nProcs);
  if (!isComplex<Cmpt>()) {
    label ir, ic;
    scalar dv;
#ifdef LABEL_INT64
    while (fscanf(fp, "%ld %ld %lf\n", &ir, &ic, &dv) != EOF) {
#else
    while (fscanf(fp, "%d %d %lf\n", &ir, &ic, &dv) != EOF) {
#endif
      Cmpt v = static_cast<Cmpt>(dv);
      label r = static_cast<label>(ir), c = static_cast<label>(ic);

      if (r >= myBeginInd + !zeroBased && r < myEndInd + !zeroBased) {
        Atuple.push_back(Triplet<Cmpt>(r, c, v));
      }
      if (r != c) {
        switch (s) {
          case SKEWSYMMETRIC:
            if (c >= myBeginInd + !zeroBased && c < myEndInd + !zeroBased) {
              Atuple.push_back(Triplet<Cmpt>(c, r, -v));
            }
            break;
          case SYMMETRIC:
            if (c >= myBeginInd + !zeroBased && c < myEndInd + !zeroBased) {
              Atuple.push_back(Triplet<Cmpt>(c, r, v));
            }
            break;
          case HERMITIAN:
            if (c >= myBeginInd + !zeroBased && c < myEndInd + !zeroBased) {
              Atuple.push_back(Triplet<Cmpt>(c, r, v));
            }
            break;
          default:
            break;
        }
      }
    }
  } else {
    scalar vr = 0, vi = 0;
    label ir, ic;
#ifdef LABEL_INT64
    while (fscanf(fp, "%ld %ld %lf %lf\n", &ir, &ic, &vr, &vi) != EOF) {
#else
    while (fscanf(fp, "%d %d %lf %lf\n", &ir, &ic, &vr, &vi) != EOF) {
#endif
      Cmpt v = getScalar<Cmpt>(vr, vi);
      label r = static_cast<label>(ir), c = static_cast<label>(ic);
      if (r >= myBeginInd + !zeroBased && r < myEndInd + !zeroBased) {
        Atuple.push_back(Triplet<Cmpt>(r, c, v));
      }
      if (r != c) {
        switch (s) {
          case SKEWSYMMETRIC:
            if (c >= myBeginInd + !zeroBased && c < myEndInd + !zeroBased) {
              Atuple.push_back(Triplet<Cmpt>(c, r, -v));
            }
            break;
          case SYMMETRIC:
            if (c >= myBeginInd + !zeroBased && c < myEndInd + !zeroBased) {
              Atuple.push_back(Triplet<Cmpt>(c, r, v));
            }
            break;
          case HERMITIAN:
            if (c >= myBeginInd + !zeroBased && c < myEndInd + !zeroBased) {
              Atuple.push_back(Triplet<Cmpt>(c, r, myConj(v)));
            }
            break;
          default:
            break;
        }
      }
    }
  }  // namespace UNAP
  fclose(fp);

  A.lRows() = myEndInd - myBeginInd;
  A.lCols() = A.lRows();
  A.lnnz() = Atuple.size();

  if (!zeroBased)
    for (label i = 0; i < A.lnnz(); ++i) {
      (Atuple[i].r_)--;
      (Atuple[i].c_)--;
    }

  A.comm()->allReduce("sumLocalnnz", &A.lnnz(), &A.nnz(), 1, COMM_LABEL,
                      COMM_SUM);

  A.beginRow() = myBeginInd;
  A.beginCol() = myBeginInd;

  A.vtxdist().resize(nProcs + 1);
  for (label i = 0; i < nProcs; ++i) {
    if (i < remainder) {
      A.vtxdist(i) = i * (avgColNum + 1);
    } else {
      A.vtxdist(i) = i * avgColNum + remainder;
    }
  }
  A.vtxdist(nProcs) = A.tsize();

  std::stable_sort(Atuple.begin(), Atuple.end(), PATCH::compareTri<Cmpt>);
  A.ptr() = new label[A.lRows() + 1];
  A.ind() = new label[A.lnnz()];
  A.val() = new Cmpt[A.lnnz()];
  label row = -1, count = 0;
  for (label i = 0; i < A.lnnz(); i++) {
    A.val(i) = Atuple[i].v_;
    A.ind(i) = Atuple[i].c_;
    label newRow = Atuple[i].r_;
    if (newRow != row) {
      A.ptr(count) = i;
      row = newRow;
      count++;
    }
  }
  A.ptr(A.lRows()) = A.lnnz();

  A.comm()->finishTask("sumLocalnnz");
  A.splitOffDiag();
  A.setupSpMVBuffers();
  A.diag();
}

template <typename Cmpt>
CSRMatrix<Cmpt>::~CSRMatrix() {
  DELETE_POINTER(diagPtr_);
  DELETE_POINTERS(valTopo_);

#ifdef UNATPP
  DELETE_POINTER(topo_);
#endif
}

template struct csr<scalar>;
template struct csr<char>;
template struct csr<std::complex<scalar> >;
template class CSRMatrix<scalar>;
template class CSRMatrix<char>;
template class CSRMatrix<std::complex<scalar> >;
template CSRMatrix<scalar>* product(const CSRMatrix<scalar>& A,
                                    const CSRMatrix<scalar>& B);
template CSRMatrix<char>* product(const CSRMatrix<char>& A,
                                  const CSRMatrix<char>& B);
template CSRMatrix<std::complex<scalar> >* product(
    const CSRMatrix<std::complex<scalar> >& A,
    const CSRMatrix<std::complex<scalar> >& B);
template void readMatrixMarket<scalar>(CSRMatrix<scalar>& A,
                                       const char* filename, bool zeroBased);
template void readMatrixMarket<std::complex<scalar> >(
    CSRMatrix<std::complex<scalar> >& A, const char* filename, bool zeroBased);

void CSRGraph::permuteLocal(const std::vector<label>& order,
                            const std::vector<label>& iorder, label low,
                            label high) {
  std::vector<label> ptr(ptr_.size()), ind(ind_.size());
  label nnz = 0;
  label n = size();
  for (label i = 0; i < n; i++) {
    label lb = ptr_[iorder[i]];
    label ub = ptr_[iorder[i] + 1];
    for (label j = lb; j < ub; j++) {
      label c = ind_[j];
      ind[nnz++] = (c >= low && c < high) ? order[c - low] : c;
    }
    ptr[i + 1] = nnz;
  }
  std::swap(ptr_, ptr);
  std::swap(ind_, ind);
}

void CSRGraph::sortRows() {
  label n = size();
  for (label r = 0; r < n; r++) {
    std::stable_sort(&ind_[ptr_[r]], &ind_[ptr_[r + 1]]);
  }
}

}  // namespace UNAP
